#include "vrt_controller.h"

VRTController::VRTController(QStringList options, QObject *parent) :
    QObject(parent),
    m_options(options)
{

    qDebug() << "**************************************************";
    double solarIntensity, zenith, azimuth;
    solarIntensity = zenith = azimuth = 0;
    double wl = 0.555;

    QString savePath;

    bool isRayleighAtmosphere = false;

    double tau = -1.0;
    double albedo = 0;

    int numLayer = 0;

    double aerosolStdHeight = 3000; //m
    double aerosolLowerLimit = 0; //m
    double aerosolUpperLimit = 10000; //m

    QString aerosolScatTemplate;

    struct complex m1 = {1.5, 0};
    struct complex m2 = {1.33, 0};
    double rMean1 = 0.2;
    double sigma1 = 0.45;
    double rMean2 = 0.5;
    double sigma2 = 0.45;
    double numRatio1_2 = 15.0;

    bool isScatMatrixOnly = false;

    if(m_options.indexOf("-wl") != -1){
        int i = m_options.indexOf("-wl");

        bool getNumOk = false;
        wl = m_options.at(i+1).toDouble(&getNumOk);
        if(!getNumOk){
            qDebug() << "Error wavelength parameter: " + m_options.at(i+1);
            exit(-1);
        }
        qDebug() << "Input wavelength:" << wl;
    }

    if(m_options.indexOf("-caculate_scat_matrix_only") != -1){
        isScatMatrixOnly = true;
    }

    if(m_options.indexOf("-solar") != -1){
        int i = m_options.indexOf("-solar");

        bool getNumOk = false;
        solarIntensity = m_options.at(i+1).toDouble(&getNumOk);
        zenith = m_options.at(i+2).toDouble(&getNumOk);
        azimuth = m_options.at(i+3).toDouble(&getNumOk);
        zenith = zenith * M_PI / 180;
        azimuth = azimuth * M_PI /180;
        if(m_options.indexOf("-Rayleigh_atmosphere") == -1){
            //Aerosol atmosphere
            int iz = qRound((zenith-ANGLE_STEP/2)/ ANGLE_STEP);
            int ia = qRound(azimuth / ANGLE_STEP);
            zenith = ANGLE_STEP/2 + iz * ANGLE_STEP;
            azimuth = ia * ANGLE_STEP;
        }

        if(!getNumOk || solarIntensity <= 0){
            qDebug() << "Error solar parameter: " + m_options.at(i+1) + m_options.at(i+2) + m_options.at(i+3);
            exit(-1);
        }
        qDebug() << "Input solar parameters(angle rounded):" << solarIntensity << zenith*180/M_PI << azimuth*180/M_PI;
    } else{
        if(!isScatMatrixOnly){
            qDebug() << "must specify solar parameters: -solar Intensity zenith azimuth";
            exit(-1);
        }
    }

    if(m_options.indexOf("-path") != -1){
        int i = m_options.indexOf("-path");
        savePath = m_options.at(i+1);
    }

    if(m_options.indexOf("-numlayer") != -1){
        int i = m_options.indexOf("-numlayer");
        bool getNumOk = false;

        numLayer = m_options.at(i+1).toInt(&getNumOk);
        if(!getNumOk || numLayer < 1){
            qDebug() << "Error number_of_layer parameter: " + m_options.at(i+1);
            exit(-1);
        }
        qDebug() << "Input number_of_layer parameter:" << numLayer;
    } else{
        if(!isScatMatrixOnly){
            qDebug() << "must specify number_of_layer: -numlayer VALUE";
            exit(-1);
        }
    }

    if(m_options.indexOf("-tau") != -1){
        int i = m_options.indexOf("-tau");
        bool getNumOk = false;

        tau = m_options.at(i+1).toDouble(&getNumOk);
        if(!getNumOk || tau <= 0){
            qDebug() << "Error optical depth parameter: " + m_options.at(i+1);
            exit(-1);
        }
        qDebug() << "Input optical depth parameter:" << tau;
    } else{
        if(!isScatMatrixOnly){
            qDebug() << "must specify optical depth: -tau VALUE";
            exit(-1);
        }
    }

    if(m_options.indexOf("-albedo") != -1){
        int i = m_options.indexOf("-albedo");
        albedo = m_options.at(i+1).toDouble();
        if(albedo >= 0 && albedo <= 1)
            qDebug() << "Input ground albedo parameter:" << albedo;
        else{
            qDebug() << "Input ground albedo parameter error:" << m_options.at(i+1);
        }
    }

    if(m_options.indexOf("-Rayleigh_atmosphere") != -1){
        isRayleighAtmosphere = true;
        goto PARAMETER_PARSE_OVER;
    }

    if(m_options.indexOf("-aerosol_height") != -1){
        int i = m_options.indexOf("-aerosol_height");
        aerosolStdHeight = m_options.at(i+1).toDouble();
        aerosolLowerLimit = m_options.at(i+2).toDouble();
        aerosolUpperLimit = m_options.at(i+3).toDouble();
        qDebug() << "Aerosol height profile, standard height:" << aerosolStdHeight << ", lower height limit:" << aerosolLowerLimit << ", and upper height limit:" << aerosolUpperLimit;
    }

    if(m_options.indexOf("-aerosol_scat_matrix") != -1){
        int i = m_options.indexOf("-aerosol_scat_matrix");
        aerosolScatTemplate= m_options.at(i+1);
        qDebug() << "aerosol_scat_matrix: " << aerosolScatTemplate;
    }

    if(m_options.indexOf("-radius") != -1){
        int i = m_options.indexOf("-radius");
        rMean1 = m_options.at(i+1).toDouble();
        rMean2 = m_options.at(i+2).toDouble();
        qDebug() << "Average radius of small and big particle: " << QString().sprintf("{%.3f, %.3f}", rMean1, rMean2);
    }

    if(m_options.indexOf("-sigma") != -1){
        int i = m_options.indexOf("-sigma");
        sigma1 = m_options.at(i+1).toDouble();
        sigma2 = m_options.at(i+2).toDouble();
        qDebug() << "Standard deviation of radius of small and big particle: " << QString().sprintf("{%.3f, %.3f}", sigma1, sigma2);
    }

    if(m_options.indexOf("-m1") != -1){
        int i = m_options.indexOf("-m1");
        m1.re = m_options.at(i+1).toDouble();
        m1.im = m_options.at(i+2).toDouble();
        if(m1.im > 0){
            qDebug() << "Error: imagin part of refractive index must be 0 or negative, while you input:" << m_options.at(i)+" "+m_options.at(i+1)+" "+m_options.at(i+2);
            exit(-10);
        }else
            qDebug() << "Averaged complex refrative index for small particle: " + QString().sprintf("{%.3f, %f}", m1.re, m1.im);
    }
    if(m_options.indexOf("-m2") != -1){
        int i = m_options.indexOf("-m2");
        m2.re = m_options.at(i+1).toDouble();
        m2.im = m_options.at(i+2).toDouble();
        if(m2.im > 0){
            qDebug() << "Error: imagin part of refractive index must be 0 or negative, while you input:" << m_options.at(i)+" "+m_options.at(i+1)+" "+m_options.at(i+2);
            exit(-10);
        }else
            qDebug() << "Averaged complex refrative index for big particle: " + QString().sprintf("{%.4f, %f}", m2.re, m2.im);
    }

    if(m_options.indexOf("-num_ratio_small_big") != -1){
        int i = m_options.indexOf("-num_ratio_small_big");
        numRatio1_2= m_options.at(i+1).toDouble();
        qDebug() << "Optical depth ratio of small aerosol particles to big ones(small/big): " << numRatio1_2;
    }

PARAMETER_PARSE_OVER:
    qDebug() << "***************************************************";

    AtmosLayer::setWaveLength(wl);
    if(isRayleighAtmosphere){
        RayleighTopLayer layerTop("Top layer", 300000, 400000);

        Layers midLayers;
        for(int i = 1; i < numLayer - 1; i++){
            midLayers.append(new RayleighAtmosphere(QString().sprintf("Mid layer%d",i), 20000, 300000));
        }

        RayleighGroundLayer layerGround("Ground layer", 0, 20000);

        layerGround.setGroundAlbedo(albedo);
        layerTop.setResultSaveMode(SAVE_ALL);
        layerGround.setResultSaveMode(SAVE_ALL);
        if(!savePath.isEmpty())
            layerTop.setSavePath(savePath);

        m_layers.append(&layerTop);
        m_layers.append(midLayers);
        m_layers.append(&layerGround);
        connectLayers(m_layers);

        for(int i = 0; i < m_layers.size(); i++){
            AtmosLayer *lay = m_layers[i];
            lay->setTau(tau/m_layers.size());
        }
        layerTop.RayleighCoeffInit();
        layerTop.recvSolar(solarIntensity, zenith, azimuth);
        QCoreApplication::exec();


    } else{ //******Aerosol Atmosphere
        //******Top layer is always Rayleigh scatter atmosphere.***
        RayleighTopLayer layerTop("Top layer", aerosolUpperLimit, 1000000);

        //*****Now for aerosol atmosphere***************
        //first the heights of layers from lower to upper
        QList<double> heights(AtmosLayer::expHeightProfileCaculate(aerosolStdHeight, aerosolLowerLimit, aerosolUpperLimit, numLayer));

        double (*smatrix)[3][3] = new double[N_SCAT][3][3];
        double (*pbsdfAero)[N_ZENITH][2*N_PI][3][3] = new double[N_ZENITH][N_ZENITH][2*N_PI][3][3];

        if(aerosolScatTemplate.isEmpty()){
            AerosolAtmosphere::aerosolInitTwoPeak(rMean1, sigma1, m1,
                                                  rMean2, sigma2, m2,
                                                  numRatio1_2, smatrix);
            if(m_options.indexOf("-caculate_scat_matrix_only") != -1){
                exit(0);
            }
        } else{
            AerosolAtmosphere::readAerosolScatMatrix(aerosolScatTemplate,
                                                     smatrix);
        }
        AerosolAtmosphere::averageAerosolMatrix(smatrix, pbsdfAero);
        Layers midLayers;
        for(int i = heights.size() - 1; i > 1; i--){
            midLayers.append(new AerosolAtmosphere(QString().sprintf("Mid layer%d",i), heights[i-1], heights[i], smatrix, pbsdfAero, tau/numLayer));
        }
        if(!midLayers.isEmpty())
            midLayers.at(0)->setResultSaveMode(SAVE_UP);

        AerosolGroundLayer layerGround("Ground layer", heights[0], heights[1], smatrix, pbsdfAero, tau/numLayer);

        layerGround.setGroundAlbedo(albedo);
        layerTop.setResultSaveMode(SAVE_ALL);
        layerGround.setResultSaveMode(SAVE_ALL);
        if(!savePath.isEmpty())
            layerTop.setSavePath(savePath);

        m_layers.append(&layerTop);
        if(!midLayers.isEmpty())
            m_layers.append(midLayers);
        m_layers.append(&layerGround);
        connectLayers(m_layers);
        connect(&layerGround, SIGNAL(caculationOver()),
                this, SLOT(overActions()));

        layerTop.RayleighCoeffInit();
        layerTop.recvSolar(solarIntensity, zenith, azimuth);

        m_solarIntensity = solarIntensity;
        m_solarZen = zenith;
        m_solarAzi = azimuth;

        m_ground = &layerGround;

        QCoreApplication::exec();
    }
}

void VRTController::overActions()
{
}

//Connect a series of layers up-down signal/slots.
//Assumming layer[0] is the toppest layer.
void VRTController::connectLayers(Layers l)
{
    int i = 0;
    AtmosLayer *upper, *lower;

    while(i < l.length() - 1){
        upper = l[i];
        lower = l[i+1];
//qDebug() << "upper lower's name:" << upper->getLayerName() << lower->getLayerName();

        QObject::connect(upper, SIGNAL(emitLightDown(QByteArray)),
                         lower, SLOT(recvScatInDown(QByteArray)));
        QObject::connect(upper,
                         SIGNAL(emitSolarLight(double,double,double)),
                         lower,
                         SLOT(recvSolar(double,double,double)));

        QObject::connect(lower, SIGNAL(emitLightUp(QByteArray)),
                         upper, SLOT(recvScatInUp(QByteArray)));
        i++;
    }
}

AtmosLayer *VRTController::getTopLayer()
{
    return m_layers.first();
}

double VRTController::getSolarIntensity()
{
    return m_solarIntensity;
}

double VRTController::getSolarZenith()
{
    return m_solarZen;
}

double VRTController::getSolarAzimuth()
{
    return m_solarAzi;
}
